Individual tree segmentation of airborne and UAV LiDAR point clouds based on the watershed and optimized connection center evolution clustering

Abstract Light detection and ranging (LiDAR) data can provide 3D structural information of objects and are ideal for extracting individual tree parameters, and individual tree segmentation (ITS) is a vital step for this purpose. Various ITS methods have been emerging from airborne LiDAR scanning (ALS) or unmanned aerial vehicle LiDAR scanning (ULS) data. Here, we propose a new individual tree segmentation method, which couples the classical and efficient watershed algorithm (WS) and the newly developed connection center evolution (CCE) clustering algorithm in pattern recognition. The CCE is first used in ITS and comprehensively optimized by considering tree structure and point cloud characteristics. Firstly, the amount of data is greatly reduced by mean shift voxelization. Then, the optimal clustering scale is automatically determined by the shapes in the projection of three different directions. We select five forest plots in Saihanba, China and 14 public plots in Alpine region, Europe with ULS or ALS point cloud densities from 11 to 3295 pts/m2. Eleven ITS methods were used for comparison. The accuracy of tree top detection and tree height extraction is estimated by five and two metrics, respectively. The results show that the matching rate (R match) of tree tops is up to 0.92, the coefficient of determination (R 2) of tree height estimation is up to .94, and the minimum root mean square error (RMSE) is 0.6 m. Our method outperforms the other methods especially in the broadleaf forests plot on slopes, where the five evaluation metrics for tree top detection outperformed the other algorithms by at least 11% on average. Our ITS method is both robust and efficient and has the potential to be used especially in coniferous forests to extract the structural parameters of individual trees for forest management, carbon stock estimation, and habitat mapping.


| INTRODUC TI ON
Forests, as a vital part of terrestrial ecosystems, play an important role in global climate change and biodiversity Seidl et al., 2017). It is challenging to conduct resource surveys of forests, especially at the individual tree scale. In the past, forest resource surveys often relied on field measurements, which were time-consuming and laborious. In recent years, remote sensing data have been increasingly applied to forestry. 2D optical images have been used to estimate forest morphological parameters (e.g., canopy cover and leaf area index) (Korhonen et al., 2017). However, these data are unable to retrieval three dimensional (3D) structural information of trees (Zheng et al., 2021). Light detection and ranging (LiDAR) data provide 3D structural information of objects and are ideal for extracting individual tree parameters of forests (Lefsky et al., 2002). There are two main categories of LiDAR for extracting individual tree parameters: ground-based and air-based. Groundbased LiDAR, such as terrestrial LiDAR scanning (TLS), has a high distance accuracy of the measurement and denser points within the limited extent, which is suitable for delicate structural parameter extraction at the plot scale (Burt et al., 2019;Tao et al., 2015). Air-based LiDAR including airborne LiDAR scanning (ALS) and unmanned aerial vehicle LiDAR scanning (ULS) can be applied to survey 3D information in a bigger region than TLS with a little lower points density.
Considering that ALS and ULS can acquire the 3D structural characteristics of trees on a large scale in complex terrain conditions, they are often used in forest survey (Guo et al., 2020).
These methods are usually more efficient, but the part of the information will inevitably be lost when the 3D point clouds are converted into 2D rasters (Zhen et al., 2016). In addition, CHMs or DSMs may also have pits, which dramatically affect the accuracy of the segmentation algorithm (Yang et al., 2019;Zhang et al., 2020).
The point-based methods directly utilize primitive or voxelized point clouds for ITS, such as point cloud region growing Lu et al., 2014), layer stacking (Ayrey et al., 2017), k-means (Lindberg et al., 2014), and graph cut (Lindberg et al., 2014;Williams et al., 2019). These methods can better use the 3D structure information of the point cloud data and further improve segmentation accuracy (Zhen et al., 2016). However, these methods also suffer from complex parameters, poor generalizability, or low efficiency.
The joint methods combine the first two in the hope of achieving a better result. For example, Tochon et al. (2015) combined the watershed and k-means algorithms to ITS in conifer and broadleaf forests. Reitberger et al. (2009) first extracted the trunk using the watershed algorithm and then used the extracted trunk as a priori knowledge of normalized cut. The joint methods combine the advantages of the first two categories of methods and therefore can improve the segmentation accuracy, but will also inherit both the disadvantage of the raster-and point cloud-based methods. In some studies, data from ALS and ULS have not been distinguished because of the similarity of their data collection principles (Yun et al., 2021). But in fact they differ significantly in point density. The point density of ALS is typically limited to 10 points/m 2 , while the point density of ULS can range from 10 to t1000 points/m 2 depending on the flight altitude and sensor characteristics (Kellner et al., 2019;Lu et al., 2014). As a result, ULS usually contains more detailed information than ALS. ITS studies for ULS have been conducting to achieve better segmentation result. For example, Wallace et al. (2014), Balsi et al. (2018) and Yin and Wang (2019) used ULS for ITS in homogenous forest. Jaskierniak et al. (2021) develop a bottom-up approach of ITS for mixed species eucalypt forests. Although these studies have get good results, the forest scenes are homogenous or specific.
Several critical issues about the presented ITS methods of ALS and ULS are summarized as follows: (1) There is an urgent need to propose more general and flexible methods that are not specific to data sources or forest types. Vauhkonen et al. (2012) compared six different ITS methods and found that the forest structure strongly affected the performance of all algorithms. Wang et al. (2016) found that point density was a highly influential factor in the performance of the methods that use point cloud data. Robust methods that are not sensitive to point density (both suit for ALS and ULS) and can be applied to coniferous, broadleaf, and mixed forests are rarely seen in the current studies. (2) There is an urgent need to propose methods that are specific to certain challenging forest types or scenarios. Dense vegetation, undulating terrain, differences in canopy shape and size, etc. can make it difficult to ITS. It is necessary to analyze the mechanism of the impact of special scenarios on ITS and propose targeted solutions. For example, the issue of omission (under-segmentation) is a big challenge for most ITS methods for dense forests (Table 1). A summary about under-and over-segmentation percentages of some ITS methods

T A X O N O M Y C L A S S I F I C A T I O N
Landscape ecology is listed in Table 1. According to the study of Li et al. (2012), when the tree stem density increases from 0.05 to 0.07 trees/m 2 , the percentage of omission greatly increases from 15% to 29% even in conifer forests. Broadleaf and mixed forests even have bigger omission fractions than conifers because of the complex structures and various species of trees. The reason for these results is that there is a severe mutual shading effect among the trees in the dense forest. Therefore, methods that make full use of the detailed information in the point cloud are needed. Joint ITS methods take the advantages of both the high efficiency of the raster-based methods and the high accuracy of the point-based methods, which have better development prospects.
The basic idea of the joint methods is to use the raster-based methods for initial segmentation and then the point-based methods for fine segmentation. Many point clustering algorithms in pattern recognition can be used for fine segmentation, such as k-means (Lindberg et al., 2014), mean-shift (Dai et al., 2018), and graphbased algorithms (Lindberg et al., 2014;Williams et al., 2019).
However, these algorithms directly rely on the input parameters, and different parameters may yield very different results (Geng & Tang, 2020). Therefore, it is necessary to develop a robust clustering algorithm that does not depend excessively on the input parameters.
In this study, we propose a new joint individual tree segmentation algorithm coupled with the watershed and optimized connection center evolution algorithm. Firstly, we use a pit-free canopy height model to implement initial segmentation based on the watershed (WS) algorithm, which has the advantages of high efficiency. Secondly, we introduce a new clustering algorithm called connection center evolution (CCE), which extends the concept of the number of paths in graph theory to the case of arbitrary real numbers and can automatically skip the unreasonable number of clusters (Geng & Tang, 2020). and then fine segmentation based on the optimized CCE algorithm, which reduced data amount by voxelization and determines the optimal clustering scale by different planar projections.
The motivation of this study is to provide individual tree attributes such as height and location for the construction of largescale digital forestry. Therefore, a general and efficient ITS method is expected. For this purpose, ALS and ULS data from different forest types, such as coniferous, broadleaf and mixed forests, with different point cloud densities were used and validated by location and tree height. This paper is organized according to the following structure. In Section 1, we introduce the overview of our study site and datasets and describe how the data are preprocessed.
The basic principle and framework of our method are explained in Section 2. In Section 3, the results and analysis are displayed.
The discussion and conclusion are explained in Sections 4 and 5, respectively.

| Plots
We selected five forest plots for the validation (Figure 1). P1 is a deciduous broadleaf forest plot (birch); P2 is a mixed forest plot containing deciduous coniferous and evergreen coniferous and deciduous broadleaf (mixed with aspen, larch, Mongolian pine, spruce, and birch); P3 is a deciduous coniferous forest plot (larch); P4 and P5 are both evergreen coniferous forests (including spruce and Mongolian pine, respectively). The area of the plots is 30 m × 30 m or 50 m × 50 m, and the average tree density of all 5 plots is 0.10 trees/ m 2 . The specific information of these plots is shown in Table 2.

| Data acquisition and preprocessing
The data in this study include both point cloud data and field measurement data. Point cloud data in each plot were obtained by ULS and TLS devices. ULS data were used to test the ITS methods, while the combination of TLS and field measurement data is used to obtain accurate reference locations for each tree. It is extremely difficult to measure the height of a large number of single trees in the field, especially for our study area where the tree height is usually greater than 10 m. Therefore, we merged the ULS and TLS data and then manually extracted the tree height of each tree as reference. The reference tree height and location were also used to evaluate the ITS methods.

| Acquisition of LiDAR and field data
The specific data include the following three types.  Table 2. (2) TLS point clouds: A Riegl VZ-1000 terrestrial laser scanner was used to obtain multistation scanning data at the sampling center and corners in order to relieve the occlusion issue. Depending on plot size and canopy characteristics, 9-17 scanning stations were set up.
(3) Field data: Fieldwork was also carried out in August 2022. We used HI TARGET Qstar 8 Mobile GPS to locate the center points of the plots. The location of each tree in all plots was checked and corrected by manual field surveys according to tree locations extracted from TLS (see Section 2.2.2). We did not use GPS to locate each tree because of the large uncertainty in positioning in the understory.

| Data preprocessing
Data preprocessing for the ULS, TLS, and field trunk position data includes the four-step operations, which are illustrated intuitively in Figure 2. (1)

| Benchmark airborne LiDAR point clouds
To demonstrate the applicability of the ITS method in different study areas, forest types, and point cloud densities, we used a benchmark airborne LiDAR point cloud dataset with individual tree inventory data in the Alpine Space, Europe (Eysn et al., 2015).
This dataset includes 14 different plots located in four European countries and can be downloaded from the NEWFOR website (https://www.newfor.net/). The detailed descriptions of these plots are shown in Table 3. Due to the low point cloud densities in these plots, the resolutions of CHMs and DTMs generated by ALS point clouds were set to 0.2 m.

| The method
Our ITS method consists of three main components: pit-free CHM generation, initial segmentation using the WS, and fine segmentation using the optimized CCE. The implementation is shown in Figure 3.

| Pit-free CHM generation
The pit-free CHM can eliminate the pits and thus reduce oversegmentation (Yang et al., 2019

| Initial segmentation using the WS
The WS is an image region segmentation method, which takes the similarity with the neighboring pixels as an essential reference in the segmentation process so that the pixels with similar spatial locations and similar grayscale values (height value in the CHM) are connected to form a closed contour (Wang et al., 2004). Here, the lidR tools developed by Roussel et al. (2020) are used to implement the WS and get the initial segmentation results. There are two input parameters: height tolerance (denoted as tolerance) and neighborhood search radius (denoted as ext). Tolerance represents the minimum height of the object in the units of image intensity between its highest point (seed) and the point where it contacts another object (checked for every contact pixel). If the height is smaller than the tolerance, the object will be combined with one of its neighbors, which is the highest. Ext represents the radius of the neighborhood in pixels for the detection of neighboring objects. A higher ext value smoothes out small objects. Figure 4 shows an example of the ITS results by the WS. Consistent with Table 1, the method suffers from significant under-segmentation.

| Fine segmentation using the optimized CCE
In this part, the CCE algorithm is optimized and used for fine segmentation. The CCE first constructs the similarity matrix between each point, then performs the power multiplication operation on the similarity matrix continuously, and finally determines the aggregation center and the number of clusters by comparing the element sizes of the similarity matrix after each power operation.
The CCE is considered an efficient and elegant clustering algorithm in Pattern Recognition (Geng & Tang, 2020  2. Similarity matrix construction. For each "individual tree" point clouds that has been initially segmented and voxelized, we construct the point-to-point distance matrix D: the variable related to the distance between point p i and point p j . Vr is the vertical distance correction factor (value range is 0-1), which is introduced to consider the incompleteness of the ULS/ALS point clouds in the lower part of the tree canopy due to occlusion. n i and n j are the weights of the two voxels, that are used to maintain the consistency of the voxel space with the original point clouds.
Next, the distance matrix (D) can be converted to the similarity matrix (S) by the Gaussian kernel function (Geng & Tang, 2020) as follows: where σ is an empirical coefficient that controls the size of the Gaussian kernel function. The element s i,j represents the similarity between p i and p j .
The similarity matrix is similar in concept to the adjacency matrix, but the elements of the similarity matrix can be real numbers.
Typically, the elements themselves are the most similar, so the diagonal elements of the similarity matrix are maximal.
(1) Finally, the similarity matrix needs to be normalized as follows: where D is the degree matrix of S and d i is the degree of the ith point (p i ).
3. CCE clustering. First, the power operation is performed on the normalized similarity matrix to obtain the following k-order connectivity: The entry (s k i,j ) of the kth power (S k ) of the similarity matrix (S) is defined as the k-order connectivity between p i and p j (denoted as con (k) p i , p j ). In particular, the diagonal entry s k i,i is defined as the korder connectivity of point p i (denoted as con (k) p i , p i ). For each k, the k-order relative connectivity of all points can be calculated, and the clustering centers will be determined according to the following rules: If one point satisfies Equation (6), it will be a connection center of the graph and is defined as a k-order clustering center of the data.
After the clustering centers are determined, the relative connectivity (rcon (k) (i, j)) is calculated according to Equation (7), and the clustering rules (p * ) are determined according to Equation (8). If we have m clustering centers p c i c i ∈ {1, 2, … , n} and i = 1, 2, … , m , for any point p j , it will be assigned to p*, where p* satisfies Equation (8).
For some datasets, for different values of k, we may obtain the same clustered data but with slightly different clustering results. In this situation, we can retain the optimal clustering results by introducing the normalized cut as follows: where P l represents the complement of P l in P and Vol P l is the sum of k-order connectivity between all points in P l and all points in P.

4.
Automatic determination of optimal scale. According to the CCE clustering, the clustering situation of different scales can be determined. When k = 1, each point is a clustering center, which is the most microscopic case. As the value of k increases, more points will be grouped together, which is the macroscopic case. We need to determine that the clustering results of the optimal scale and correctly segment individual trees. For this purpose, we project each scale of clustering result point clouds to the X-Y, X-Z and Y-Z plane, respectively ( Figure 5). Then, we determine whether the following three inequalities hold in each of the three projection planes: F I G U R E 5 Diagram of segmented point clouds projected to X-Y, X-Z, and Y-Z planes. The blue filled graph shows the approximate outline of the tree point cloud projected onto the X-Y and Y-Z planes. Crown X and Crown Y are the crown widths along the X-axis and Y-axis direction. X Z max is the x value of the point with maximum z projected onto the X-Z plane. Y Z max is the corresponding parameter on the Y-Z plane. Three constraint inequalities are list on Corresponding projection planes.
Crown X and Crown Y are the crown widths along the X-axis and Y-axis direction. X Z max is the x value of the point with maximum z projected onto the X-Z plane, X min and X max are the maximum and minimum x values in all points projected to this plane, respectively. Y Z max , Y min , and Y max are the corresponding parameters on the Y-Z plane. Finally, we filter the clustering results that satisfy the above conditions. If there are multiple candidate results, the one with the maximum number of the candidates will be selected as the best.
Equations (10)-(12) ensures that the segmented tree shape is rational. Equation (10) requires that the larger of the crown width in the X and Y directions does not exceed three times that of the smaller, and canopies that exceed this limit are rare in nature. The sensitivity of the parameters of Equations (11) and (12) is analyzed in Section 4.1.
The input parameters of our method are summarized in Table 4. In addition, for plantations with trees of relatively similar growth, we followed the postprocessing method proposed by Pang et al. (2021). If the distance between two adjacent individual trees is less than the average crown diameter of the corresponding plot and the elevation of these two trees is less than 10 m, they will be merged into an individual tree. The average crown width is calculated using the segmented point clouds by our ITS algorithm. The watershed algorithm and data processing also involve the corresponding parameters. The sensitivity analysis of these parameters is not addressed in this study, as it has been previously analyzed by corresponding studies (Pang et al., 2016;Wang et al., 2004).
Finally, each tree height and location are automatically extracted by calculating the height and geographical coordinates of the highest point.

| Accuracy assessment
LiDAR point clouds with tree labels are output after applying the ITS method. Then, horizontal location and tree height are matching to the field reference data. The matching method started from the highest detected tree and searched for the reference trees that satisfied the height and distance criterion as match candidates. If a farther candidate showed a better height difference, then it became a better match. This process was repeated until all detected trees have been checked. If the closest one with the smallest height difference is the matched detection tree previously, these two trees will be treated as a matched pair (Pang et al., 2021). The matching criterion is described by Eysn et al. (2015). Eventually, a series of matching parameters are calculated. TP (true positive) is the number of correctly segmented trees; FN (false negative) is the number of trees not segmented but assigned to a nearby tree (omission error or under-segmentation); FP (false-positive) is the number of trees that did not exist but were segmented from the point cloud (commission error or over-segmentation).
We select extraction rate (R extraction ), matching rate (R match ), commission rate (R commission ), omission rate (R omission ), and F score (F) as evaluation metrics (Eysn et al., 2015;Li et al., 2012). Here are the expressions.
R match , R extraction , and F are the main assessment metrics, related to the overall accuracy, and the closer they are to 1, the higher the accuracy of the ITS algorithm. R omission and R commission are secondary assessment metrics to measure the degree of under-and over-segmentation, and the closer they are to 0, the less under-or over-segmentation. The above metrics are used for tree top detection, and tree height estimation is evaluated by coefficient of determination (R 2 ) and root mean square error (RMSE).

Vr
The vertical distance correction factor Reducing the influence of incompleteness of the ULS/ALS point clouds in the lower part of the tree canopy due to occlusion.

Empirical coefficient related to the Gaussian kernel
Controlling the size of the Gaussian kernel function when converting the distance matrix into the similarity matrix.

TA B L E 4
Description of two input parameters in our ITS method.
However, for P1, it is difficult to evaluate the segmentation results because the tree tops of the broadleaf cannot be observed clearly on one hand, on the other hand, the point cloud density in this plot is relatively low.
The quantitative assessment results are presented in Table 5.
Overall, the segmentation accuracy is fine with an average match rate and F-score greater than 0.7. However, there is some oversegmentation, especially in P3 and P5 with a relatively lower R omission . We checked ULS and TLS data carefully and found that the conifers in Sahanba, especially the larch, are prone to trunk bifurcation. A case is shown in Figure 7 to illustrate the phenomenon of trunk bifurcation. The phenomenon can be clearly seen in the TLS point clouds (Figure 7b,c). However, the tree trunk is not clearly visible through the ULS point clouds due to the occlusion issue, which causes it to look similar to two trees (Figure 7a).

| Tree height accuracy evaluation
The accuracy of tree height extraction is evaluated by comparing the reference with the matched tree heights. As seen in Figure 8, all the results are well except for P1. P3 and the 14 plots of benchmark dataset are the best with R 2 = .94, although the RMSE of the benchmark dataset is 1.667 m. The results for P2, P4, and P5 are relatively well, with R 2 = .79 (.74 for P2) and RMSE < 1 m. In general, our method can accurately extract the tree height of coniferous and mixed forests. For broadleaf forests, especially on slopes, the precise extraction of tree height requires more effort.

| Comparison with existing methods
To evaluate our approach more comprehensively, we choose three classical ITS methods for comparison, including the WS (Wang et al., 2004), mark-controlled watershed (denoted as MCWS) (Chen et al., 2006), and point cloud region growing segmentation (denoted as PCS) . The WS and PCS are implemented through the lidR tool (Roussel et al., 2020), and the MCWS implemented through Digital-Forestry-Toolbox (https://mpark an.github. io/Digit al-Fores try-Toolb ox/). Due to the high densities of the ULS point clouds in P2-P5, the PCS cannot be executed effectively.

F I G U R E 6
Diagram of the tree top detection of five plots in Saihanba. The black discrete points are the reference tree top coordinates, and the normalized point clouds with different colors are the ITS results.

TA B L E 5
Results of treetop detection of P1-P5 using our ITS algorithm.  Table 7. Compared with the other three methods, our method gives the best results. Table 8 and Figure 9 show the ITS results of 14 public plots in the benchmark dataset. Compared with the WS, MCWS, and PCS, our method gives the best matching rate. Although the F-score by our algorithm is 0.02 lower than that by MCWS, our matching rate is 0.16 higher. Our method also gives the best matching rate compared to methods #1-#8 described by Eysn et al. (2015). Of these methods, WS, #5, #6, and ours matched more than 50%. All the four methods give over-segmentation result, while ours is at the medium level.

| Sensitivity analysis and parameter settings
For P1, the result is relatively poor with the R 2 = .5. There are three reasons for this: (1) the average slope of this plot is 30°, so the point cloud normalization will cause distortion of the trees (Khosravipour et al., 2015). (2) there is distortion of the trunk of birch due to the natural environment; (3) there is no obvious top of broadleaf trees, which is different from coniferous trees. So it is difficult to accurately detect tree tops even visually. The above factors cause errors in both field measurements and algorithm estimation. Figure 10 shows the TLS point cloud data of P1 and clearly confirms the three analyses above.
P1 was the most complex plot in this study, with complex topographic conditions, the highest tree stem density, irregular canopy shape, and relatively low point cloud density. Therefore, it was used for the sensitivity and parameter settings analysis. For the optimized CCE, the optimal clustering scale is determined by Equations (10)-(12). With Equation (10), it is ensured that the shape of the canopy is reasonable and unreasonably flattened canopy is removed. With Equations (11) and (12), the distance between the top and edge is determined by projection in two directions, and then, the minimum distance threshold is set to ensure that the top is located near the center of the canopy. Table 9 demonstrates the effect of the minimum distance threshold setting on the results in P1. If no minimum distance is set (or a small value, e.g., 1/16 crown diameter), oversegmentation will be very serious. However, if this threshold is set too large (e.g., 1/4 crown diameter), many trees will not be segregated, especially for broadleaf forests with inconspicuous tree tops.
Therefore, this threshold was set to 1/8 crown diameter to ensure its applicability in both coniferous and broadleaf forests.
There are two input parameters in our algorithm. The vertical distance correction factor, Vr, is to be considered for ULS/ALS point cloud clustering. In our study, Vr is set to 1/6 according to the best results given by Pang et al. (2021). The empirical coefficient related F I G U R E 9 Tree top detection results of 14 public plots in Europe using our and other eight methods. #1-#8 correspond to the methods described by Eysn et al. (2015).
to the Gaussian kernel, σ, is was analyzed in our study. The variation of extraction rate, matching rate, commission rate, omission rate, and F score with σ 2 is shown in Figure 11. It can be seen that these five evaluation metrics are very stable, indicating that our algorithm is robust.
To fairly compare various ITS methods, the same canopy structure related parameters were set in all test plots (Table 10). These parameters are either program default parameters or determined by reference to previous studies. For the parameter related to the point cloud density, that is, the resolution of the CHM, we set this parameter to 0.2 m for ALS generation and 0.1 m for ULS generation. For the MCWS, the relationship between tree height and canopy radius is required. However, field measurements are difficult to obtain sufficient accuracy and enough data, so we refer to the formulas by Popescu and Wynne (2004) (See Table 10). In the previous section, the MCWS gave poor results in many plots. This is due to the inappropriate relationship between the tree height and crown radius within the plots, and not the algorithm itself.
The properties of different types of trees should be complex, but due to field measurements constraints, only three fixed formulas are given for broadleaf, coniferous, and mixed forests (Popescu & Wynne, 2004).

| Efficient implementation
With the development of LiDAR hardware technology, high quality and density ULS/ALS LiDAR point clouds are emerging. Therefore, ITS algorithms are also expected to be able to process data efficiently. Thanks to the initial segmentation using the watershed and the mean shift voxelization, the execution speed of the CCE algorithm has been greatly improved. The processing speed of the improved CCE algorithm was tested with the configuration of a Core Intel(R) Core(TM) i7-8700 CPU@3.20GHz Processor, 40 GB RAM, an NVIDIA GeForce GTX 1660 graphics card and the Microsoft Windows 10 operating system. The results are shown in Table 11.
We did not compare the original CCE algorithm because there

TA B L E 9
Setting of the minimum distance threshold from the top to the edge for the trees in P1.

F I G U R E 11
Variations of the tree top detection assessment metrics with σ 2 .
F I G U R E 1 0 Side view of the TLS point cloud data of P1.
was not enough memory in the device we used for the method to run successfully. For the P3 with 50 × 50 m with a density of over 1500 pts/m 2 , the time to run the algorithm is within 6 min. Our ITS method has the potential to meet the upcoming era of massive point clouds.

| CON CLUS ION
Individual tree segmentation using ALS or ULS data is still a challenge due to the complexity of forest structure. In this paper, we proposed a new individual tree segmentation method, which consists of the WS algorithm, and the optimized CCE algorithm. We optimized the CCE algorithm to make it more efficient, and the optimal segmentation scale can be determined automatically by taking into account the structural characteristics of the canopy. The new ITS method can take full advantages of the efficient of the WS and the accuracy of CCE algorithm. Additionally, the new method is robust for the complex plots and insensitive for the parameters. Tree coordinates and heights are extracted and output directly automatically.
Validation at five different forest types of plots in China and 14 public plots in Europe showed the accuracy of both treetop detection and tree height estimation. Compared with the other 11 individual tree segmentation methods, our method gives better results.
Through sensitivity analysis for input parameters, we find that the algorithm is robust. Efficient processing speed enables it to meet the high-density point clouds of 4000 pts/m 2 . Our method is both practical and applicable and can be used to extract the structural parameters of individual trees over large areas for forest management, carbon stock estimation, and habitat mapping.

TA B L E 11
Program runtime in different plots.

CO N FLI C T O F I NTE R E S T S TATE M E NT
We declare no conflicts of interest with this research.

DATA AVA I L A B I L I T Y S TAT E M E N T
The source code can be downloaded freely from https://github.com/ liyi-rs/ITS_WS_ICCE. Benchmark dataset can be available from the NEWFOR website (https://www.newfor.net/).